module force
use omplib
use technical
use general
use mod_2dflu
implicit none
save
private
double precision,allocatable,dimension(:) :: fx,fy
character (len=labellen) :: ftype='gravity-x',fprofile='none',fgravx_profile='mid-channel-zero-mean'
double precision :: famp=0.02,fchannel_bottomy=2.5,fchannel_jump=8.
integer :: fdirn=0
double precision :: kf=2.,fphase=0.
public :: read_force_pars,initialize_force,add_rspace_force,add_cont_force_fourier
namelist /force_pars/ &
  ftype,famp,fprofile,fgravx_profile,kf,fchannel_bottomy,fchannel_jump,fphase
contains
!**********************************************************!
subroutine read_force_pars(unit,ierr)
  integer, intent(in) :: unit
  integer, intent(inout), optional :: ierr
!----------
  if (present(ierr)) then
    read(unit,NML=force_pars,IOSTAT=ierr)
  else
    read(unit,NML=force_pars)
  endif
  rewind(unit)
endsubroutine read_force_pars
!-------------------------------------------------------------
subroutine initialize_force()
  select case(ftype)
    case('zero')
      call force_zero
    case('kolmogorov-x')
      call force_kolmogorov_x
    case('kolmogorov-y')
      call force_kolmogorov_y
    case('kolmogov-fourier')
        lfourier_forcing=.true.
        lrspace_forcing=.false.
   case('gravity-x')
     lfourier_forcing=.false.
     lrspace_forcing=.true.
      call force_gravity_x
    case default 
      call fatal_error('initialize_force','If you select force module, select a force')
  endselect
  call write_force()
endsubroutine initialize_force
!-------------------------------------------------------------
subroutine force_zero()
  integer :: ix
  write(*,*) ' No external force'
  write(*,*) 'setting lfourier_forcing to false'
  lfourier_forcing=.false.
  write(*,*) 'and lrspace_forcing to false'
  lrspace_forcing=.false.
endsubroutine force_zero
!-------------------------------------------------------------
subroutine force_kolmogorov_x()
  integer :: ix
  write(*,*) ' Kolmogorov forcing, x dependent'
  write(*,*) 'setting lfourier_forcing to false'
  lfourier_forcing=.false.
  write(*,*) 'and lrspace_forcing to true'
  lrspace_forcing=.true.
  write(*,*) 'f_y = famp*sin(kf*x+fphase)'
  fdirn=1
  call allocate_force()
  fx=0.
  do ix=1,n1
    fy(ix) = famp*sin(kf*xx(ix)+fphase)
  enddo
endsubroutine force_kolmogorov_x
!-------------------------------------------------------------
subroutine force_kolmogorov_y()
  integer :: iy
  write(*,*) ' Kolmogorov forcing, y dependent'
  write(*,*) 'setting lfourier_forcing to false'
  lfourier_forcing=.false.
  write(*,*) 'and lrspace_forcing to true'
  lrspace_forcing=.true.
  write(*,*) 'f_x = famp*sin(kf*y+fphase)'
  fdirn=2
  call allocate_force()
  do iy=1,n2
    fx(iy) = famp*sin(kf*yy(iy)+fphase)
  enddo
  fy=0.
endsubroutine force_kolmogorov_y
!-------------------------------------------------------------
subroutine force_gravity_x()
  integer :: ix,iy
  double precision :: ycord,ymid,yabs
  write(*,*) 'Gravity acting along x direction'
  write(*,*) 'But may also depend on y'
  select case(fgravx_profile)
  case('none')
   fdirn=0
   call allocate_force()
   fy=0.
   fx = famp  
  case('mid-channel-zero-mean')
!    we assume that the channel is symmetric about the
!    middle of the box.
!    fx=famp*0.5*(1 - Tanh(8*(Abs(y) - fchannel_bottomy)) - <fx>
   fdirn=2
    call allocate_force()
    fy=0.
    ymid=(yy0+ymax)/2
    do iy=1,n2
      ycord=yy(iy)
      yabs=abs(ycord-ymid)
      fx(iy) = famp*0.5*(1 - tanh(fchannel_jump*(yabs - fchannel_bottomy) )  )
    enddo
    fx = fx -sum(fx)*dx/length
  case default
    call fatal_error('force_gravity_y','you must select a profile')
  endselect
endsubroutine force_gravity_x
!-------------------------------------------------------------
subroutine allocate_force()
  select case(fdirn)
  case(0)
    allocate(fx(1))
    allocate(fy(1))
  case(1)
    allocate(fx(n1))
    allocate(fy(n1))
  case(2)
    allocate(fx(n2))
    allocate(fy(n2))
  case default
    call fatal_error('allocate_force','fdirn not selected')
  endselect
endsubroutine allocate_force
!-------------------------------------------------------------
subroutine write_force()
  integer :: iwrite
  open(unit=11,file='force.ascii',status='unknown')
  select case(fdirn)
  case(0)
    write(11,*) fx,fy
  case(1)
    do iwrite=1,n1
      write(11,*)xx(iwrite),fx(iwrite),fy(iwrite)
    enddo
  case(2)
    do iwrite=1,n2
      write(11,*)yy(iwrite),fx(iwrite),fy(iwrite)
    enddo
  case default
    call fatal_error('write_force','fdirn not recognized')  
  endselect
  close(11)
endsubroutine write_force
!-------------------------------------------------------------
subroutine add_rspace_force(wXu,ix,iy)
  integer, intent(in) :: ix,iy
  double precision,dimension(2) :: wXu 
!
  select case(fdirn)
  case(0)
    wXu(1) = wXu(1) + fx(1)
    wXu(2) = wXu(2) + fy(1)
  case(1) 
    wXu(1) = wXu(1) + fx(ix)
    wXu(2) = wXu(2) + fy(ix)
  case(2)
    wXu(1) = wXu(1) + fx(iy)
    wXu(2) = wXu(2) + fy(iy)
  case default
    call fatal_error('add_rspace_force','fdirn not selected')
  endselect
!  
endsubroutine add_rspace_force
!-------------------------------------------------------------
subroutine add_cont_force_fourier(k1,k2,fomegak)
  use omplib
  double precision :: force_amp
  integer :: k1,k2
  double precision :: fomegak
  force_amp=famp/scale
  if ((k1.eq.kf).and.(k2.eq.0)) then
    fomegak=force_amp*kf/2.0d0
!    call get_thread_id(id)
!    write(*,*) 'I am thread',id,fomegak,k1,k2,force_amp,kforce
  else
    fomegak=0.0d0
  endif
!
endsubroutine add_cont_force_fourier
!**********************************************************!
endmodule force
